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ABSTRACT 



'Living fossils', a phrase first coined by Darwin, are defined as species with limited 
recent diversification and high morphological stasis over long periods of evolu- 
tionary time. Morphological stasis, however, can potentially lead to diversification 
rates being underestimated. Notostraca, or tadpole shrimps, is an ancient, globally 
distributed order of branchiopod crustaceans regarded as 'living fossils' because their 
rich fossil record dates back to the early Devonian and their morphology is highly 
conserved. Recent phylogenetic reconstructions have shown a strong biogeographic 
signal, suggesting diversification due to continental breakup, and widespread cryptic 
speciation. However, morphological conservatism makes it difficult to place fossil 
taxa in a phylogenetic context. Here we reveal for the first time the timing and tempo 
of tadpole shrimp diversification by inferring a robust multilocus phylogeny of Bran - 
chiopoda and applying Bayesian divergence dating techniques using reliable fossil 
calibrations external to Notostraca. Our results suggest at least two bouts of global 
radiation in Notostraca, one of them recent, so questioning the validity of the 'living 
fossils' concept in groups where cryptic speciation is widespread. 
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INTRODUCTION 

There has been much debate about the tempo and mode of the diversification of life 
{Eldredge & Gould, 1972; Reznick & Ricklefs, 2009; Rhodes, 1983). Recently, this debate 
has been informed by dating using relaxed molecular clocks and diversification analyses; 
techniques which have revealed disparate patterns of speciation with early bursts (Burbrink 
& Pyron, 2010), recent radiations (Nagalingum et al, 2011) and density dependency 
(Phillimore & Price, 2008) being demonstrated. One extreme and often controversial 
pattern of diversification is found in 'living fossils', a concept introduced by Charles 
Darwin in On the Origin of Species when dealing with the perplexing nature of the platypus 
and lungnsh, relicts of once diverse groups (Darwin, 1859). Since Darwin's first use, the 
'living fossil' term has been applied to groups which appear to have diversified little and are 
morphologically stable over long periods of evolutionary time, with examples including 
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cycads, tuatara, coelacanths, horseshoe crabs and Ginkgo biloba. However, morphological 
stasis can obscure the patterns of species diversification, and recent time-calibrated 
phylogenetic analyses of some 'living fossils' have indeed revealed that extant species are in 
fact only recently diverged (Kano et al, 2012; Nagalingum et al, 201 1 ) . 

Notostraca, or tadpole shrimps, is an ancient, globally distributed order of branchiopod 
crustaceans with a rich fossil record dating back to the early Devonian (Fayers & Trewin, 
2002). The order has two extant genera, Triops and Lepidurus, in the family Triopsidae, 
with a yet undefined number of species. The nomenclature and systematic position 
of some ancient extinct Notostraca lineages, is, however, problematic (Hegna & Dong, 
2010). This is partly because tadpole shrimps have maintained an extremely conserved, 
yet complex, bauplan with extant species indistinguishable from fossils of Triops from the 
Triassic {Gall & Grauvogel-Stamm, 2005; Gore, 1986; Trusheim, 1938) and of Lepidurus in 
the Jurassic (Barnard, 1929; Houghton, 1924). This striking morphological conservatism 
has led them to be referred to as 'living fossils' (Fryer, 1988; King & Hanner, 1998; 
Mantovani et al, 2008). 

Phylogenetic reconstructions of extant Notostraca show a strong biogeographic signal 
(Mathers et al, 2013; Vanschoenwinkel et al, 2012). In Triops, species complexes are largely 
restricted to single continents, while Lepidurus lineages show high levels of endemism 
(Rogers, 2001), patterns that suggest ancient radiation with diversification through 
continental break-up. However, the extreme morphological conservatism of this order 
hampers both the taxonomy of extant species and the phylogenetic placement of fossil 
taxa, with little known about the timing and tempo of notostracan diversification. Genetic 
analyses have revealed widespread cryptic species (King & Hanner, 1998; Korn et al, 2010; 
Korn & Hundsdoerfer, 2006; Macdonald, Sallenave & Cowley, 2011; Vanschoenwinkel etal, 
2012), further illustrating the difficulty of inferring past and present diversity. To address 
this difficulty we infer a robust phylogeny of all known notostracan species from both 
extant genera and seven branchiopod outgroups. Our analysis uses all available Notostraca 
sequence data for seven genes, and Bayesian relaxed clock dating techniques, with multiple 
branchiopod fossil calibrations, to estimate divergence times. 

MATERIALS AND METHODS 

Species delimitation 

As Notostraca is known to contain cryptic species complexes (e.g. King & Hanner, 1998), 
and in order to follow the same criterion for species selection for the multilocus analysis, 
we delimited species using a generalised mixed Yule coalescent (GMYC) model (Pons et al, 
2006) fitted to an ultrametric phylogeny based on all available cytochrome oxidase I (COI) 
sequences from GenBank. The 270 sequences were aligned with Muscle (Edgar, 2004) and 
phylogeny estimated with BEAST vl .7.4 (Drummond et al, 2012) under a constant popu- 
lation size coalescent tree model and GTR +F substitution model. A strict molecular clock 
was used with the substitution rate fixed to 1 to provide branch lengths relative to an arbi- 
trary time scale. The MCMC chain was run for 9,000,000 iterations with the first 500,000 
iterations removed as burnin. Effective sample sizes (ESS's) of parameters (all greater 
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than 200) and appropriate burnin were checked using Tracer vl.5 (Rambaut & Drummond, 
2007). From this a maximum clade credibility tree using median heights was made. 
We then fitted a single threshold GMYC model to the COI tree to delimit species from 
populations. A total of 34 species of Notostraca were identified in this analysis (Table SI). 

Multilocus phylogenetic analysis of Branchiopoda 

We constructed a multilocus alignment containing representatives of all known species 
of Notostraca. Single representatives of each phylogenetic species identified by the 
GMYC analysis were selected for inclusion in our phylogenetic analysis. In addition, 
four Notostraca lineages (T. gadensis, T. cf. granarius Tunisia, L. bilobatus andL. cryptus), 
which did not have COI data available, but were represented by other genes, were also 
used in our multilocus phylogenetic analysis. The species status of these lineages has been 
confirmed in regional studies of cryptic diversity (King & Hanner, 1998; Korn et ah, 201 0; 
Korn & Hundsdoerfer, 2006; Rogers, 2001). We also included seven representatives of the 
other branchiopod orders so that robust fossil calibrations could be applied for the dating 
analysis. We included sequences for the genes 12S, 16S, 28S, cytochrome oxidase I (COI), 
Elongation Factor 1 -alpha (EF1), RNA Polymerase II and Glycogen Synthase (see Table S2 
for Accession Numbers). 

Sequences were aligned using MUSCLE (Edgar, 2004) with final adjustments by 
eye. Introns in the nuclear protein coding genes were identified and removed based 
on alignment with available Notostraca mRNAs. Translation was checked in MEGA 5 
(Tamura et ah, 2011). Overall, sequences for 45 taxa (38 notostracan and 7 branchiopod 
outgroups) were concatenated for analysis with the alignment containing 5793 positions 
and 52% missing data (Table S2; the alignment file is available in Dryad DOI 10.5061/ 
dryad.77bt2). 

Optimum partitioning schemes and substitution models for our phylogenetic analysis 
and divergence time estimation were identified using PartitionFinder (Lanfear et ah, 2012). 
PartitionFinder uses a heuristic search algorithm, starting with a fully partitioned analysis 
(gene and codon position where appropriate), and identifies the best fit partitioning 
scheme and substitution models based on Bayesian Information Criterion (BIC). Due to 
the restricted model choice available in RAxML (Stamatakis, 2006) we conducted separate 
PartitionFinder analyses for the phylogenetic analysis and divergence time estimation. 
For the phylogenetic analysis we restricted model choice options to GTR or GTR +r, 
whereas for divergence time estimation we allowed models to be selected from the full 
suite available in BEAST. We excluded models with proportion of invariant sites (+1) 
as rate heterogeneity is accounted for by the gamma shape parameter (+T). Optimum 
partitioning schemes and substitution models for both analyses are given in Tables S3 
and S4. 

Branchiopod phylogeny was estimated using Bayesian and maximum likelihood (ML) 
methods with partitions and substitution models set to those identified by PartitionFinder 
(Table S3). Bayesian analysis was performed with MrBayes v3.2 (Ronquist et ah, 2012). 
Model parameters between partitions were unlinked. Two independent MCMC chains 
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Table 1 Fossils used to calibrate divergence time analysis in BEAST. Age constraints are treated 
numbers indicate phylogenetic placement of fossil calibrations in Fig. 1. 


as hard bounds unless otherwise stated. Node 


Node 


Fossil taxa 


Geological period 


Minimum age 

(Mya) 


Maximum age 

(Mya) 


Reference 


1 


Oldest Bilateria eg. Kimberella 


Ediacaran 




558 


Fedonkin, Simonetta & 
Ivantsov (2007) 


1 


Rehbachiella 


Upper Cambrian 




500 (soft max) 


Walofiek (1995) 


1 


Undescribed anostracan 


Base Ordovician 


488 




Harvey, Velez & 
Butterfield (2012) 


2 


Castracollis 


Pragian, Early Devonian 


410 




Fayers & Trewin (2002) 


3 


Ebullitiocaris elatus 


Carboniferous 


300 




Womacketal. (2012) 


4 


Daphnia and Ctenodaphnia sp. 


Jurassic/Cretaceous 


145 




Kotov& Taylor (2011) 



were run for 10,000,000 iterations each, sampling every 5,000 iterations. The first 25% 
of each run was discarded as burnin with the remaining samples pooled and used to 
create a maximum clade credibility tree. Maximum likelihood phylogenetic analysis was 
performed using RAxMLHPC-PTHREADS v7.0.4 (Stamatakis, 2006). An initial ML search 
using GTR + T4 was performed onto which 100 rapid bootstraps were drawn. 

Bayesian relaxed clock divergence dating 

We estimated Bayesian divergence times with BEAST vl .7.4 (Drummond et al, 2012) using 
an uncorrelated lognormal relaxed clock (Drummond et al, 2006) and a Yule speciation 
prior. XML files for all BEAST runs were created using BEAUTi vl.7.4 (Drummond 
et al, 2012). Topology was constrained to that of the unconstrained RAxML analysis. 
We used the best fit partitioning scheme identified by PartitionFinder (Table S4) and 
estimated substitution model parameters independently for each partition. Initial runs 
were conducted using substitution models identified by PartitionFinder, however, this 
resulted in poor mixing of some GTR model parameters for partitions 1, 2 and 5, so 
subsequent runs were performed using a simpler HKY +F model for these partitions. 

Five branchiopod fossils representing the oldest known occurrences of their respective 
crown groups were used to calibrate the molecular clock with minimum age constraints 
(Table 1). Lognormal prior distributions were used to specify the level of uncertainty in 
the placement of these fossil calibrations as they reflect the likely scenario that the true date 
of divergence of a given node was some time before the earliest known fossil belonging 
to that clade (Ho & Phillips, 2009). Lognormal distributions have three parameters - 
mean, standard deviation and offset. We set the offset to correspond to the minimum 
age of the node as determined by the fossil record, we then specified mean and standard 
deviations that resulted in 95% of the distribution falling between the age of the fossil 
and the age of the next oldest fossil (at a lower taxonomic level) for that group. This gives 
a prior distribution, which assigns the majority of the probability close to the age of the 
oldest known fossil, but gives a long tail to account for uncertainty in the proximity of 
the fossil to the true date of divergence. As Bayesian divergence dating benefits from at 
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least one maximum age constraint we conservatively constrained the root of the tree to a 
maximum age of 558 Mya, the age of the oldest bilaterian fossil (Fedonkin, Simonetta & 
Ivantsov,2007). 

We ran two independent BEAST MCMC chains for 50,000,000 iterations, sampling 
every 5000 iterations. Ten million iterations were removed as burnin from each run. 
Convergence of the two runs and the ESS of parameter estimates (all greater than 
250) where assessed using Tracer vl.5 (http://tree.bio.ed.ac.uk/software/tracer/). A 
posterior sample of 8000 trees from one of the runs was used to construct a maximum 
clade credibility time tree for Notostraca and our selected outgroups. The XML file 
used to run the BEAST divergence time analysis can be downloaded from Dryad 
DOI 10.5061/dryad.77bt2. 

Diversification analysis 

Patterns of diversification through time within Notostraca were investigated using LASER 
(Rabosky, 2006) based on the BEAST time tree with outgroups pruned. Using LASER we 
compared constant rate and variable rate speciation models using likelihood ratio tests 

(Table S5). 

RESULTS AND DISCUSSION 

The most complete taxon sampling to date coupled with the inclusion of multiple nuclear 
and mitochondrial markers allowed us to generate a robust phylogeny of 38 extant 
Notostraca species. ML and Bayesian inference gave highly congruent phylogenetic trees 
with most branches highly supported (Figs. SI and S2). The recovered relationships 
between branchiopod orders are in agreement with recently published arthropod 
phylogenies (Regier et al, 2010; von Reumont et al., 2012), providing a solid platform 
for divergence dating analysis. 

Our robust time-calibrated phylogeny of Branchiopoda (Fig. 1) clearly shows that 
Notostraca has a pattern of diversity incompatible with Darwin's original usage of the term 
'living fossil' as relics of once diverse groups, and importantly reveals cryptic patterns of 
diversification. Our analysis - using outgroup fossil calibrations - estimates an ancient 
divergence of Triops and Lepidurus during the Jurassic, 184 Mya (95% confidence interval 
132-259 Mya), which agrees with the earliest fossils assigned to Lepidurus (Barnard, 1929; 
Gand, Garric & Lapeyrie, 1997), and with a sister relationship of Notostraca to the extinct 
order Kazacharthra of the Late Triassic/Early Jurassic (Olesen, 2009). This initial radiation 
of extant Notostraca was not, however, due to continental break up as the timing and 
pattern of diversification within the genera substantially postdates the break-up of Pangaea 
160-138 Mya (Scotese, 2001). Furthermore, the current species distributions of Triops and 
Lepidurus are likely to have resulted from a second global radiation of the order, possibly 
following considerable levels of extinction. We conclude this because fossil Notostraca, 
attributed to Triops and Lepidurus, have been found in modern day North and South 
America, Europe, Africa and Antarctica, implying a global distribution by the early Jurassic 
(Gall & Grauvogel-Stamm, 2005; Gand, Garric & Lapeyrie, 1997; Garrouste, Nel & Gand, 
2009; Gore, 1986; Haughton, 1924; Trusheim, 1938). Yet, our LASER analysis shows a 
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■ Palearctic 


■ Afrotropic 


■ Nearctic 


■ Australasia 


■ Indomalaya 






T. cf. australiensis (Paynes Find) 
T. cf. australiensis sp. 2 
T. sp. B {Paroo) 

T. cf.australiensis(W alga Rock) 

T. cf. australiensis(Gibb Rock) 

T. sp. A (Paroo) 

T. cf. australiensis sp. 3 

T. cf. australiensis (Baladonia Rock) 

T. cf australiensis (Ayers Rock) 

T. cf. australiensis sp. I 

T. newberryi 

T. cf. longicaudatus sp. 2 

T. cf. longicaudatus sp. I 

T. cf australiensis (Lake Carey) 

T. cf. granarius (South Africa) 

T. cf. granarius (Tunisia) 

T. cf. granarius (Namibia) 

T. cf. granarius (Russia) 

T. cf. granarius (Japan) 

T. vicentinus 

T. emeritensis 

T. mauritanicus mauhtanicus 
T. simplex 

T. cf. mauritanicus (East Spain) 

'/'. gadensis 

T. haeticus 

T. cancriformis 

L. cf couesii (Apulia) 

L. cf. couesii (Canada) 

L. couesii (Sardinia) 

L. arcticus 

L. viridis 

L. cryptus 

L. lemmoni 

L. bilobatus 

L. packardi 

L. opus 

L. lubbocki 

Eulimnadia sp. 

Limnadia lenticularis 

Daphnia pulex 

Daphnia magna 

Lynceus sp. 

Anemia sp. 

Streptocephalus sp. 



Camb. 


Ordovician 


Silurian 


Devonian 


Carboniferous 


Permian 


Triassic 


Jurassic 


Cretaceous 


Paleogene 


Neog. 


Paleozoic 


Mesozoic 


Cenozoic 



0 Mya 



Figure 1 Time calibrated phylogeny of 38 notostracan species and seven branchiopod outgroups. Numbers at nodes correspond to the fossil 
calibrations given in Table 1. Nodes with black circles have ML Bootstrap support values greater than 70 and posterior probabilities greater than 95 
from the RAxML and MrBayes analyses respectively. Error bars show the 95% confidence intervals of divergence times. Colour coded squares show 
the known geographic distribution of each species. 



significant increase in the rate of diversification of Notostraca about 73 Mya (Fig. 2), 
close to the time of the K-Pg mass extinction event. It is this second radiation that resulted 
in the current global distribution of extant Triops and Lepidurus. 

The almost synchronous radiation of Triops and Lepidurus (Fig. 1) suggests that a 
common factor may have triggered diversification of the two genera. The diversification 
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Figure 2 Diversification of Notostraca through time. Arrows indicate the timing and direction of 
shifts in rate of diversification inferred by LASER. N is the number of species. The best fit ML model 
of diversification identified three distinct rates of diversification during the evolutionary history of 
Notostraca with an increase in speciation rate 73 Mya followed by a decrease 6 Mya. 



of modern birds - widely involved in dispersal in aquatic invertebrates (Green et al, 2005; 
van Leeuwen et al, 2012) — coincided with the initiation of the notostracan radiation 
(Pacheco et al, 201 1), and may have facilitated the long distance dispersal and subsequent 
diversification of Notostraca. Indeed, the geographical distribution of extant taxa (Fig. 1) 
suggests several instances of intercontinental dispersal. For example, the colonisation 
of North America from Australia could have resulted from dispersal events during bird 
migration. Such long distance dispersal and colonisation events might also have been 
facilitated by the flexible nature of sexual systems found within Notostraca (Mathers 
et al, 2013). Indeed, the evolution of androdioecy - a sexual system where males and 
hermaphrodites coexist (Weeks, 2012; Zierold, Hanfling & Gomez, 2007; Zierold et al, 
2009) - from gonochorism appears to have favoured postglacial recolonisation in the 
species Triops cancriformis (Zierold, Hanfling & Gomez, 2007) and in Notostraca as a whole 
(Mathers et al, 2013). 

The concept of 'living fossils' has been a controversial one as it has often been 
interpreted to imply a lack of evolutionary change, even against evidence of molecular 
evolutionary change (Avise, Nelson & Sugita, 1994; Casane & Laurenti, 2013). Our 
divergence dating analysis has shown that tadpole shrimps can be regarded as 'living 
fossils' only on the grounds of morphological conservatism, not on the basis of limited 
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diversification or relict status. Instead, throughout their long evolutionary history, 
notostracans have undergone multiple global radiations and high species turnover. Recent, 
time calibrated, phylogenetic analysis of other traditional 'living fossils' such as cycads 
(Nagalingum et al, 2011), nautiloids (Wray et al, 1995), horseshoe crabs (Obst et al, 
2012) and monoplacophorans (Kano et al, 2012), have also revealed that extant species are 
more recently diverged than suggested by fossil data alone. We therefore caution against 
drawing conclusions about patterns of diversification based on fossil data alone in groups 
where widespread morphological conservatism may obscure rampant cryptic speciation. 
Furthermore, our results help clarify the term 'living fossils', putting important questions 
into focus. Namely, is such morphological conservatism, in the face of evolutionarily 
recent diversification and radiation, best accounted for by unchanging selection or by 
developmental genetic constraints? 
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